﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace 高斯正反算
{
    public class Guass2
    {
        protected int zone_wide = 3;
        protected double semi_major = 6378137.0;           // major axis
        protected double semi_minor = 6356752.3142451793;           // minor axis
        protected double scale_factor = 1.0;    // scale factor
        protected double central_meridian;      //Center longitude (projection center)
        protected double lat_origin = 0;        //center latitude
        protected double e0, e1, e2, e3, e, es, esp;    // eccentricity constants
        protected double ml0;                   //small value m
        protected double false_northing = 0.0;  //y offset in meters
        protected double false_easting;         //x offset in meters

        #region ICoordinate Members
        public int ZoneWide
        {
            get { return zone_wide; }
            set { zone_wide = value; }
        }

        public void DD2DMS(double DecimalDegree, out int Degree, out int Minute, out double Second)
        {
            Degree = (int)DecimalDegree;
            Minute = (int)((DecimalDegree - Degree) * 60.0);
            Second = Math.Round((DecimalDegree * 60 - Degree * 60 - Minute) * 60 * 100) / 100.0;

        }

        //public double DistanceOfTwoPoints(double lng1, double lat1, double lng2, double lat2, GaussSphere gs)
        //{
        //    double radLat1 = Degrees2Radians(lat1);
        //    double radLat2 = Degrees2Radians(lat2);
        //    double a = radLat1 - radLat2;
        //    double b = Degrees2Radians(lng1) - Degrees2Radians(lng2);
        //    double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) +
        //    Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2)));
        //    s = s * (gs == GaussSphere.WGS84 ? 6378137.0 : (gs == GaussSphere.Xian80 ? 6378140.0 : 6378245.0));
        //    s = Math.Round(s * 10000) / 10000;
        //    return s;
        //}

        protected static double Degrees2Radians(double deg)
        {
            return (Math.PI / 180 * deg);

        }

        protected static double Radians2Degrees(double rad)
        {
            return (180 / Math.PI * rad);
        }

        public void DMS2DD(int Degree, int Minute, double Second, out double DecimalDegree)
        {
            DecimalDegree = Degree + Minute / 60.0 + Second / 60.0 / 60.0;
        }

        public void CaleZoneIndex(double longitude, out int zoneIndex, out double centralMeridian)
        {
            if (this.zone_wide == 6)
            {
                zoneIndex = (int)(longitude / (double)this.zone_wide) + 1;
                centralMeridian = 6 * zoneIndex - 3;
            }
            else
            {
                zoneIndex = (int)((longitude - 1.5) / (double)this.zone_wide) + 1;
                centralMeridian = 3 * zoneIndex;
            }
        }

        public void GaussPrjCalculate(double longitude, double latitude, out double X, out double Y)
        {
            double lon = Degrees2Radians(longitude);
            double lat = Degrees2Radians(latitude);

            double delta_lon = 0.0;         //Delta longitude (Given longitude - center
            double sin_phi, cos_phi;        //sin and cos value
            double al, als;                 //temporary values
            double c, t, tq;                //temporary values
            double con, n, ml;              //cone constant, small m

            es = 1.0 - Math.Pow(semi_minor / semi_major, 2);
            e = Math.Sqrt(es);
            e0 = (1.0 - 0.25 * es * (1.0 + es / 16.0 * (3.0 + 1.25 * es)));
            e1 = (0.375 * es * (1.0 + 0.25 * es * (1.0 + 0.46875 * es)));
            e2 = (0.05859375 * es * es * (1.0 + 0.75 * es));
            e3 = (es * es * es * (35.0 / 3072.0));
            ml0 = semi_major * (e0 * lat_origin - e1 * Math.Sin(2.0 * lat_origin) + e2 * Math.Sin(4.0 * lat_origin) - e3 * Math.Sin(6.0 * lat_origin));
            esp = es / (1.0 - es);

            int projNo = 0;
            if (this.zone_wide == 6)
            {
                projNo = (int)(longitude / (double)this.zone_wide) + 1;
                this.central_meridian = Degrees2Radians((double)(projNo * this.zone_wide - (this.zone_wide >> 1)));
            }
            else
            {
                projNo = (int)((longitude - 1.5) / (double)this.zone_wide) + 1;
                this.central_meridian = Degrees2Radians((double)projNo * (double)this.zone_wide);
            }
            false_easting = 1000000L * projNo + 500000L;

            //int projNo = (int)(longitude / zone_wide);
            //central_meridian = Degrees2Radians(projNo * zone_wide + zone_wide / 2);
            //false_easting = 1000000L * (projNo + 1) + 500000L;


            delta_lon = Adjust_lon(lon - central_meridian);
            sin_phi = Math.Sin(lat);
            cos_phi = Math.Cos(lat);

            al = cos_phi * delta_lon;
            als = Math.Pow(al, 2);
            c = esp * Math.Pow(cos_phi, 2);
            tq = Math.Tan(lat);
            t = Math.Pow(tq, 2);
            con = 1.0 - es * Math.Pow(sin_phi, 2);
            n = semi_major / Math.Sqrt(con);
            ml = semi_major * (e0 * lat - e1 * Math.Sin(2.0 * lat) + e2 * Math.Sin(4.0 * lat) - e3 * Math.Sin(6.0 * lat));


            X = scale_factor * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 *
                (5.0 - 18.0 * t + Math.Pow(t, 2) + 72.0 * c - 58.0 * esp))) + false_easting;
            Y = scale_factor * (ml - ml0 + n * tq * (als * (0.5 + als / 24.0 *
                (5.0 - t + 9.0 * c + 4.0 * Math.Pow(c, 2) + als / 30.0 * (61.0 - 58.0 * t
                + Math.Pow(t, 2) + 600.0 * c - 330.0 * esp))))) + false_northing;
        }

        public void GaussPrjInvCalculate(double X, double Y, out double longitude, out double latitude)
        {
            double con, phi;                    //temporary angles
            double delta_phi;                   //difference between longitudes
            long i;                             //counter variable
            double sin_phi, cos_phi, tan_phi;   //sin cos and tangent values
            double c, cs, t, ts, n, r, d, ds;   //temporary variables
            long max_iter = 6;                  //maximun number of iterations

            es = 1.0 - Math.Pow(semi_minor / semi_major, 2);
            e = Math.Sqrt(es);
            e0 = (1.0 - 0.25 * es * (1.0 + es / 16.0 * (3.0 + 1.25 * es)));
            e1 = (0.375 * es * (1.0 + 0.25 * es * (1.0 + 0.46875 * es)));
            e2 = (0.05859375 * es * es * (1.0 + 0.75 * es));
            e3 = (es * es * es * (35.0 / 3072.0));
            ml0 = semi_major * (e0 * lat_origin - e1 * Math.Sin(2.0 * lat_origin) + e2 * Math.Sin(4.0 * lat_origin) - e3 * Math.Sin(6.0 * lat_origin));
            esp = es / (1.0 - es);

            int projNo = (int)(X / 1000000L);
            if (this.zone_wide == 6)
            {              
                this.central_meridian = Degrees2Radians((double)(projNo * this.zone_wide - (this.zone_wide >> 1)));
            }
            else
            {               
                this.central_meridian = Degrees2Radians((double)projNo * (double)this.zone_wide);
            }
            false_easting = 1000000L * projNo + 500000L;

            double x = X - false_easting;
            double y = Y - false_northing;

            con = (ml0 + y / scale_factor) / semi_major;
            phi = con;
            for (i = 0; ; i++)
            {
                delta_phi = ((con + e1 * Math.Sin(2.0 * phi) - e2 * Math.Sin(4.0 * phi) + e3 * Math.Sin(6.0 * phi))
                    / e0) - phi;
                phi += delta_phi;
                if (Math.Abs(delta_phi) <= 1.0e-10) break;
                if (i >= max_iter)
                    throw new ApplicationException("Latitude failed to converge");
            }
            if (Math.Abs(phi) < Math.PI / 2.0)
            {
                sin_phi = Math.Sin(phi);
                cos_phi = Math.Cos(phi);
                tan_phi = Math.Tan(phi);
                c = esp * Math.Pow(cos_phi, 2);
                cs = Math.Pow(c, 2);
                t = Math.Pow(tan_phi, 2);
                ts = Math.Pow(t, 2);
                con = 1.0 - es * Math.Pow(sin_phi, 2);
                n = semi_major / Math.Sqrt(con);
                r = n * (1.0 - es) / con;
                d = x / (n * scale_factor);
                ds = Math.Pow(d, 2);

                double lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t +
                    10.0 * c - 4.0 * cs - 9.0 * esp - ds / 30.0 * (61.0 + 90.0 * t +
                    298.0 * c + 45.0 * ts - 252.0 * esp - 3.0 * cs)));
                double lon = Adjust_lon(central_meridian + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t +
                    c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * esp +
                    24.0 * ts))) / cos_phi));
                longitude = Radians2Degrees(lon);
                latitude = Radians2Degrees(lat);
            }
            else
            {
                longitude = Radians2Degrees(Math.PI / 2 * (y < 0.0 ? -1 : 1));
                latitude = Radians2Degrees(central_meridian);
            }
        }

        private double Adjust_lon(double x)
        {
            long count = 0;
            double prjMAXLONG = 2147483647;
            double DBLLONG = 4.61168601e18;
            for (; ; )
            {
                if (Math.Abs(x) <= Math.PI)
                {
                    break;
                }
                else
                {
                    if (((long)Math.Abs(x / Math.PI)) < 2)
                        x = x - ((x < 0.0 ? -1 : 1) * Math.PI * 2.0);
                    else
                        if (((long)Math.Abs(x / Math.PI * 2.0)) < prjMAXLONG)
                        {
                            x = x - (((long)(x / Math.PI * 2.0)) * Math.PI * 2.0);
                        }
                        else
                        {
                            if (((long)Math.Abs(x / (prjMAXLONG * Math.PI * 2.0))) < prjMAXLONG)
                            {
                                x = x - (((long)(x / (prjMAXLONG * Math.PI * 2.0))) * (Math.PI * 2.0 * prjMAXLONG));
                            }
                            else
                            {
                                if (((long)Math.Abs(x / (DBLLONG * Math.PI * 2.0))) < prjMAXLONG)
                                {
                                    x = x - (((long)(x / (DBLLONG * Math.PI * 2.0))) * (Math.PI * 2.0 * DBLLONG));
                                }
                                else
                                {
                                    x = x - ((x < 0.0 ? -1 : 1) * Math.PI * 2.0);
                                }
                            }
                        }
                    count++;
                    if (count > 4) break;
                }
            }
            return (x);
        }
        #endregion
    }
}
